function bcoef = computegpdistance(X1,X2)
% Determine bhat. coefficient between multivariate time series X1 and X2

D = length(X1);
N = length(X1{1});
vec = zeros(size(X1{1},1),1);

if D == 1
    s1 = [X1{1}(:,1) X1{1}(:,2)];
    s2 = [X2{1}(:,1) X2{1}(:,2)];
    [~, vec] = gpdist(s1,s2);
    
    bcoef = median(vec(2:2:length(vec)));

else
    
    for n = 1 : N
        mu1 = zeros(1,D);
        mu2 = zeros(1,D);
        cov1 = zeros(D);
        cov2 = zeros(D);
        for d = 1 : D
            mu1(1,d) = X1{d}(n,1); mu2(1,d) = X2{d}(n,1);
            cov1(d,d) = X1{d}(n,2);
            cov2(d,d) = X2{d}(n,2);
        end
        [~,vec(n)] = bhat_gauss(mu1,cov1,mu2,cov2);
    end
    
    bcoef = mean(vec(2:2:length(vec)));
    
end